function K = getContour(I, thres)

[sy, sx] = size(I);
K = zeros(sy, sx);
P = zeros(16, 4);
P(1,:)= [0 0 0 0]; % [0 0;0 0]
P(2,:)= [0 0 0 1]; % [1 0;0 0] 
P(3,:)= [0 0 1 0]; % [0 1;0 0] 
P(4,:)= [0 0 1 1]; % [1 1;0 0] 
P(5,:)= [0 1 0 0]; % [0 0;1 0] 
P(6,:)= [0 1 0 1]; % [1 0;1 0] 
P(7,:)= [0 1 1 0]; % [0 1;1 0] ambiguous
P(8,:)= [0 1 0 1]; % [1 1;1 0] 
P(9,:)= [1 0 0 0]; % [0 0;0 1] 
P(10,:)=[1 0 0 1]; % [1 0;0 1] ambiguous
P(11,:)=[1 0 1 0]; % [0 1;0 1] 
P(12,:)=[1 0 1 0]; % [1 1;0 1] 
P(13,:)=[1 1 0 0]; % [0 0;1 1] 
P(14,:)=[0 1 0 1]; % [1 0;1 1] 
P(15,:)=[1 0 1 0]; % [0 1;1 1] 
P(16,:)=[0 0 0 0]; % [1 1;1 1]

B = I>=thres;

for y = 1:(sy+1)
	for x = 1:(sx+1)
		ul = 0;
		ur = 0;
		ll = 0;
		lr = 0;
        
		both = (x > 1) & (y > 1) & (x <= sx) & (y <= sy);
		if (x > 1 && y <= sy)
			ll = B(y, x - 1);
		end
		if (y > 1 && x <= sx)
			ur = B(y - 1, x);
		end
		if (both)
			ul = B(y - 1, x - 1);
        end
        if ((x <= sx) && (y <= sy))
            lr = B(y, x);
        end
		
		index = 8 * ul + 4 * ur + 2 * ll + lr + 1;
		contourP = P(index,:);
		
		if (contourP(1) == 1) 
			K(y - 1, x - 1) = 255; 
		end
		if (contourP(2) == 1)
			K(y - 1, x) = 255;
		end
		if (contourP(3) == 1)
			K(y, x - 1) = 255;
        end
        if (contourP(4) == 1)
            K(y, x) = (contourP(4) * 255);
        end
		
	end
end
